#!/usr/bin/env python
import sys
from astropy.io import fits
from glob import glob
import os
import numpy as np
from astropy.table import Table
if not os.path.exists("origin"):
    os.system("mkdir origin")
    os.system("mv *.fits origin")
nums = len(glob("origin/bias_outlier_15_*.fits") )
npix15 = []
npix25 = []
frac15 = []
frac25 = []
npix = []
if sys.argv[1] == '0083':
    ymin = 0
    ymax = 4616
    xmin = 0
    xmax = 9216
if sys.argv[1] == '0095':
    ymin = 4616
    ymax = 9232
    xmin = 0
    xmax = 9216
if sys.argv[1] == 'cmos':
    ymin = 0
    ymax = 9120
    xmin = 0
    xmax = 8900
for i in range(nums):
    print(i)
    bias_out_15 = fits.getdata("origin/bias_outlier_15_"+str(i+1)+".fits")[ymin:ymax,xmin:xmax]
    bias_out_25 = fits.getdata("origin/bias_outlier_25_"+str(i+1)+".fits")[ymin:ymax,xmin:xmax]
    stat_flag = fits.getdata("origin/bias_stat_flag_"+str(i+1)+".fits")[ymin:ymax,xmin:xmax]
    bias_out_15 *= stat_flag
    bias_out_25 *= stat_flag
    npix.append(np.sum(stat_flag))
    npix15.append(np.sum(bias_out_15))
    npix25.append(np.sum(bias_out_25))
    frac15.append(npix15[i]/npix[i])
    frac25.append(npix25[i]/npix[i])
    if i==0:
        count_map_15 = bias_out_15.astype('int')
        count_map_25 = bias_out_25.astype('int')
    else:
        count_map_15 += bias_out_15.astype('int')
        count_map_25 += bias_out_25.astype('int')
if not os.path.exists("remove_fix"):
    os.system("mkdir remove_fix")
    os.system("mkdir remove_fix/42")
    os.system("mkdir remove_fix/48")
    os.system("mkdir remove_fix/54")
fits.writeto("remove_fix/count_map_15.fits",data=count_map_15,overwrite=True)
fits.writeto("remove_fix/count_map_25.fits",data=count_map_25,overwrite=True)

npix15_42 = npix15.copy()
npix15_48 = npix15.copy()
npix15_54 = npix15.copy()
npix25_42 = npix25.copy()
npix25_48 = npix25.copy()
npix25_54 = npix25.copy()
frac15_42 = frac15.copy()
frac15_48 = frac15.copy()
frac15_54 = frac15.copy()
frac25_42 = frac25.copy()
frac25_48 = frac25.copy()
frac25_54 = frac25.copy()
npix15_42 = (np.array(npix15)-np.sum(count_map_15>42)).astype('int').tolist()
npix25_42 = (np.array(npix25)-np.sum(count_map_25>42)).astype('int').tolist()
npix15_48 = (np.array(npix15)-np.sum(count_map_15>48)).astype('int').tolist()
npix25_48 = (np.array(npix25)-np.sum(count_map_25>48)).astype('int').tolist()
npix15_54 = (np.array(npix15)-np.sum(count_map_15>54)).astype('int').tolist()
npix25_54 = (np.array(npix25)-np.sum(count_map_25>54)).astype('int').tolist()
frac15_42 = np.round((np.array(npix15_42)/np.array(npix)),8).tolist()
frac25_42 = np.round((np.array(npix25_42)/np.array(npix)),8).tolist()
frac15_48 = np.round((np.array(npix15_48)/np.array(npix)),8).tolist()
frac25_48 = np.round((np.array(npix25_48)/np.array(npix)),8).tolist()
frac15_54 = np.round((np.array(npix15_54)/np.array(npix)),8).tolist()
frac25_54 = np.round((np.array(npix25_54)/np.array(npix)),8).tolist()

tab_15 = Table()
tab_25 = Table()
tab_15['imgid'] = np.append(np.arange(nums).astype('int').astype('str'),np.array(['median','max']))
tab_15['npix'] = np.append(np.array(npix),np.array([-99,-99]))
tab_25['npix'] = np.append(np.array(npix),np.array([-99,-99]))
tab_15['npix15'] = np.append(np.array(npix15),np.array([np.median(npix15),np.max(npix15)]))
tab_15['frac15'] = np.append(np.array(frac15),np.array([np.median(frac15),np.max(frac15)]))
tab_25['npix25'] = np.append(np.array(npix25),np.array([np.median(npix25),np.max(npix25)]))
tab_25['frac25'] = np.append(np.array(frac25),np.array([np.median(frac25),np.max(frac25)]))
tab_15['npix15_42'] = np.append(np.array(npix15_42),np.array([np.median(npix15_42),np.max(npix15_42)]))
tab_15['frac15_42'] = np.append(np.array(frac15_42),np.array([np.median(frac15_42),np.max(frac15_42)]))
tab_25['npix25_42'] = np.append(np.array(npix25_42),np.array([np.median(npix25_42),np.max(npix25_42)]))
tab_25['frac25_42'] = np.append(np.array(frac25_42),np.array([np.median(frac25_42),np.max(frac25_42)]))
tab_15['npix15_48'] = np.append(np.array(npix15_48),np.array([np.median(npix15_48),np.max(npix15_48)]))
tab_15['frac15_48'] = np.append(np.array(frac15_48),np.array([np.median(frac15_48),np.max(frac15_48)]))
tab_25['npix25_48'] = np.append(np.array(npix25_48),np.array([np.median(npix25_48),np.max(npix25_48)]))
tab_25['frac25_48'] = np.append(np.array(frac25_48),np.array([np.median(frac25_48),np.max(frac25_48)]))
tab_15['npix15_54'] = np.append(np.array(npix15_54),np.array([np.median(npix15_54),np.max(npix15_54)]))
tab_15['frac15_54'] = np.append(np.array(frac15_54),np.array([np.median(frac15_54),np.max(frac15_54)]))
tab_25['npix25_54'] = np.append(np.array(npix25_54),np.array([np.median(npix25_54),np.max(npix25_54)]))
tab_25['frac25_54'] = np.append(np.array(frac25_54),np.array([np.median(frac25_54),np.max(frac25_54)]))
tab_15.write("result_wo_fix_15.tab",format='ipac',overwrite=True)
tab_25.write("result_wo_fix_25.tab",format='ipac',overwrite=True)

for j in range(nums):
    print(j)
    bias_out_15 = fits.getdata("origin/bias_outlier_15_"+str(j+1)+".fits")[ymin:ymax,xmin:xmax]
    bias_out_25 = fits.getdata("origin/bias_outlier_25_"+str(j+1)+".fits")[ymin:ymax,xmin:xmax]
    bias_out_15_42 = bias_out_15*(count_map_15<42).astype('uint8')
    bias_out_25_42 = bias_out_25*(count_map_25<42).astype('uint8')
    fits.writeto("remove_fix/42/bias_mask_15_42_"+str(j+1)+".fits",data=(count_map_15>42).astype('uint8'),overwrite=True)
    fits.writeto("remove_fix/42/bias_mask_25_42_"+str(j+1)+".fits",data=(count_map_25>42).astype('uint8'),overwrite=True)  
    fits.writeto("remove_fix/42/bias_outlier_15_42_"+str(j+1)+".fits",data=bias_out_15_42,overwrite=True)
    fits.writeto("remove_fix/42/bias_outlier_25_42_"+str(j+1)+".fits",data=bias_out_25_42,overwrite=True)
    bias_out_15_48 = bias_out_15*(count_map_15<48).astype('uint8')
    bias_out_25_48 = bias_out_25*(count_map_25<48).astype('uint8')
    fits.writeto("remove_fix/48/bias_mask_15_48_"+str(j+1)+".fits",data=(count_map_15>48).astype('uint8'),overwrite=True)
    fits.writeto("remove_fix/48/bias_mask_25_48_"+str(j+1)+".fits",data=(count_map_25>48).astype('uint8'),overwrite=True)
    fits.writeto("remove_fix/48/bias_outlier_15_48_"+str(j+1)+".fits",data=bias_out_15_48,overwrite=True)
    fits.writeto("remove_fix/48/bias_outlier_25_48_"+str(j+1)+".fits",data=bias_out_25_48,overwrite=True)
    bias_out_15_54 = bias_out_15*(count_map_15<54).astype('uint8')
    bias_out_25_54 = bias_out_25*(count_map_25<54).astype('uint8')
    fits.writeto("remove_fix/42/bias_mask_15_54_"+str(j+1)+".fits",data=(count_map_15>54).astype('uint8'),overwrite=True)
    fits.writeto("remove_fix/42/bias_mask_25_54_"+str(j+1)+".fits",data=(count_map_25>54).astype('uint8'),overwrite=True)
    fits.writeto("remove_fix/54/bias_outlier_15_54_"+str(j+1)+".fits",data=bias_out_15_54,overwrite=True)
    fits.writeto("remove_fix/54/bias_outlier_25_54_"+str(j+1)+".fits",data=bias_out_25_54,overwrite=True)

